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Abstract 

We consider the effects of anisotropic coherency strain (due to lattice mismatch) on phase sep- 
aration in intercalation materials, motivated by the high-rate Li-ion battery material Li a; FeP04. 
Using a phase-field model coupled to elastic stresses, we analyze spinodal decomposition (linear 
instability of the homogeneous state) as well as nonlinear evolution of the phase pattern at 
constant mean filling. We consider fully anisotropic coherency strain and focus on the case of 
simultaneous expansion and contraction along different crystal axes, as in the case of Li a; FeP04, 
which leads to tilted, striped phase boundaries in equilibrium. 

1 Introduction 

Intercalation, i.e. the reversible insertion of molecules in a host solid, is a general phenomenon 
with many applications in chemical physics. Of particular interest today is the intercalation of 
lithium ions in the active solid particles of Li-ion battery composite electrodes. Li-ion batteries 
have become essential for portable electronic devices such as cell phones and power tools and are 
advancing toward possible use in electric vehicles or renewable energy storage, due to higher energy 
density and power density than earlier battery chemistries [TH3]. In particular, Li^FePC^ has 
gained considerable attention as a leading cathode material candidate because of its high power 
density, fast charging times, long cycle life, and favorable safety characteristics [1HS]- 

Modeling ion intercalation dynamics in crystalline particles is a complex problem at the in- 
tersection of physics, chemistry, materials science, and applied mathematics [THS]. A standard 
assumption in Li-ion battery models is that diffusion is linear and isotropic within the solid [TjllO]. 
but it is becoming recognized that this is often not the case, due to the anisotropic crystalline 
nature of most solid hosts. In particular, experiments and ab initio computer simulations have 
revealed that Li diffusion in LLjFePC^ is highly anisotropic and exhibits a strong tendency for 
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phase-separation between lithiated (LiFePC^) and delithiated (FeP04) states [TTHT3] . Addition- 
ally, this FeP04/LiFeP04 phase boundary tends to reach the active (010) surface and align itself 
along {101} or {100} crystal planes [1H , I13 [ IT5]. which suggests an intimate connection between 
phase separation dynamics and the elastic properties of the crystal. In contrast, motivated by the 
first experimental paper [3j, existing battery engineering models for LizFePC^assume a spherical 
"shrinking core" of one phase being replaced by a shell of the other [16pi7j and thus do not model 
phase separation dynamics, anisotropic diffusion or coherency strain. 

Over the past five years, our group has developed a general phase- field model for intercalation ki- 
netics in nanoparticles, which is helping to unravel the complex dynamical behavior of LLjFePC^and 
interpret experimental data from microscopic first principles [18H24] , A crucial contribution of this 
work has been to consistently formulate the kinetics of electrochemical charge-transfer reactions 
at the surface of the active particles with the non-equilibrium thermodynamics of bulk ion trans- 
port and phase separation via a modified Butler- Volmer theory for concentrated solutions and 
solids [18,24 . This has led to novel predictions for intercalation dynamics in nano particles, such 
as size-dependent miscibility and spinodal gaps [21^125] . ion intercalation waves filling the crystal 
layer by layer (rather than shrinking core) at low currents [19,20,26j, and quasi-solution behavior 
and suppressed phase separation at high currents (22JI23] . 

In this work, we focus on understanding the general effects of coherency strain on bulk phase 
separation, while neglecting the coupling to electrochemical reactions. (For a complete model 
of LLjFePC^nanoparticles under applied current, see Ref. [23].) Classical phase-field models in 
materials science have included effects of elastic coherency strain for isotropic solids [34 , 35HHJH2] , 
but here we analyze phase separation with anisotropic coherency strain, focusing on the novel 
case where both contraction and expansion occur along different crystal planes, as in LLjFePC^. 
Simplified elastic models of Liz FeP04 have led to the conclusion that phase separation occurs along 
the {100} crystal planes [26.36,37], as in some experiments [15] . even though {101} phase boundaries 
are also sometimes observed [111113] . Here we show that such tilted phase boundaries result from 
simultaneous contraction and expansion in different directions. To address dynamical behavior, 
we perform a linear stability analysis on the model to predict parameter regimes for spinodal 
decomposition of the system, which provides analytical predictions for initial instabilities in the 
concentration and deformation fields. As the nonuniform states are highly nonlinear, numerical 
simulations are necessary to examine the long-term behavior of phase-separated solutions. 

2 Model Formulation 

To develop a phase-field model of this system, we begin by defining an order parameter c to be 
the dimensionless, normalized concentration of Li, where c = 1 corresponds to an entirely lithiated 
phase (LiFeP04), and c = to an entirely delithiated one (FePC^). The total free energy of the 
system described by the domain Q can be expressed as a functional of the local Li concentration 
in the Cahn-Hilliard formulation 



where /(c) is the free energy density of a homogeneous state, p is the number density of intercalation 
sites, and the gradient term (with K being a symmetric, positive definite tensor) is the next order 
term in an expansion about this state which penalizes spatial phase modulation in the system 
[271128] . For our model, we take this tensor to be constant. Li insertion into this material has been 
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observed to produce significant anisotropic strains which range from about -2% to 5% to maintain 
coherency between states [36], and hence we also take into account long-range elastic contributions 
to the free energy as in previous phase-field models [MIES]- The elastic energy is represented in 
the final term of the integrand in Equation (pQ) with the contraction between the tensors E and 
T, which are the strain and stress tensors respectively. To capture the effects of the interactions 
between lithiated and delithiated states, we take the regular solution model for the homogeneous 
energy density [27] 

/(c) = apcQ. - c) + pk B T [clog(c) + (1 - c) log(l - c)] , (2) 

where k B is the Boltzmann constant, T is the absolute temperature, and the regular solution 
interaction parameter a is the average energy density of ion intercalation. The first term in Equation 
([2]) is the enthalpic contribution which favors phase separation, while the second term is the entropic 
contribution which promotes phase mixing. The elements of the adjusted infinitesimal strain tensor 
are given by the kinematic relationship 

E = - (Vu + Vu T ) - cM, (3) 

where u is the deformation vector, and M is the lattice mismatch tensor of the lithiated and un- 
lithiated phases. This additional mismatch term is the stress- free strain when the density is c. 
Noting that the intercalation of these particles takes the delithiated orthorhombic crystal (FeP04) 
to another orthorhombic state (LiFePC^), insertion will hence only induce compressions and ex- 
pansions of the surrounding material. Consequently, the tensor M will be diagonal. Assuming the 
strains to be within the regime of linear elasticity, we can then use Hooke's law to describe the 
stress in terms of the constitutive equation 

T = C : E, (4) 

where C is the rank-4 elasticity tensor [38], which for orthorhombic crystals, can be described by 
nine distinct elements |33j . We now use this free energy to calculate the chemical potential as 



_ 1 5^ _ 1 d 

^ p 5c p dc 



f(c) + : T 



V-(KVc), (5) 



ap(l - 2c) +pk B T log (-^-), ~(E:T) = -M:T. (6) 



where 



dc " \1 — cj 2dc 

Note that as M is diagonal, only the compressional components of the stress will affect the chemical 
potential. The mass flux is then expressed as 

j = -pcBV/i, (7) 

where B is the mobility tensor. Strictly speaking, the mobility must depend on concentration, e.g. 
as B ~ (1 — c) in the classical Cahn-Hilliard regular solution model [28J, but below we will make 
simpler assumptions that have relatively little effect on the dynamics of coherency phase separation. 

Finally, the evolution of the concentration field is calculated with a conservation law within the 
cathode 

p§ + V-j = 0, (8) 
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and at any given time, a quasi-steady approximation of mechanical equilibrium is given by 

VT = 0. (9) 

Using Equations (J3J) , © and ([9]) , we can write this elastostatic condition in terms of the deformation 
field components as the system 



(10a) 



d ^ du k d (d Ul du 2 \ d f du 3 duA ^ dc 

ox f— ' oxk ay \ oy ox J oz \ ox oz J f— ' ox 

W E + Oaf + + (^1 + ^) = £ C 2t M„£, (10b) 

ay OTfc oz \ Oz dy J ox \ dy Ox J dy 

d A oufc fl /ftis 8uA d (due du 3 \_^ dc 

oxfc ax \ ax az / ay \ az ay / Oz 



fe=l v ' 1:7 v " ' fc=i 

where we are using the standard convention to reduce the elasticity tensor elements Cy&j to the 
matrix components Cy [39|. This hence gives the elastic energy contribution to the chemical 
potential 



-^-(E:T) = -M:T = V 
2 dc K ' f-f 



3 / 3u ■ \ 



(ii) 



With the bulk equations established for Li migration within the cathode, boundary conditions 
must be formulated at the electrode-electrolyte interfaces to close the model. Mass balance at the 
interface must satisfy 

n-j + fl(c,/i) = 0, (12) 

where n is the outward normal to the surface, and R(c, /i) is the potential-dependent interfacial 
reaction rate [l8 ll2~lTf2"4"] . We additionally assume the so-called variational boundary condition 

n • (KVc) = 0, (13) 

which is reasonable for our system in which concentration gradients are observed to be perpendicular 
to the particle surface [?]. Finally, we assume a uniform pressure field of strength P from the 
surrounding matrix, which gives the interfacial stress balance 

n-(T + PI) = 0. (14) 

It should be noted that we state the boundary conditions above just to show the form of the 
complete model, but we will not use them in our calculations below. As the particular focus of 
this work is a feature of the bulk dynamics, we take the infinite domain Q = M 3 , and hence only 
boundedness in the far-field of the concentration and deformation fields is necessary for the purposes 
of our analysis. Detailed studies of current dependence and finite size effects for ion intercalation 
in nanoparticles can be found in 
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3 Isotropic Coherency Strain 



To begin to understand the leading order effects of lithiation induced strains in Li x FeP04, we start 
by taking the host material to be an isotropic body. Under this assumption, the coefficients of the 
elasticity tensor reduce to 



Cn — Q 



22 



C: 



33 



1 



l-2u 

Cl2 = C21 = C23 = C32 = C13 = C31 



G, C44 — C55 

2/2 



l-2v 



G, 



G. 



(15) 
(16) 



where G is the shear modulus, and v is the Poisson ratio. We further take the coherency strains to 
be isotropic, which simplifies the mismatch tensor to the form M = M; so I, where M- lso is a scalar. 
We can then write the elastostatic condition ft9h in terms of the deformation vector as 



(1 - 2i/)V 2 u + V(V • u) = 2M iso (l + i/)Vc. 
Note that this also simplifies the elastic contribution to the chemical potential to 



ld_ 

2dc 



(E : T) = -M iso tr{T} = 2M iso G 



l + i/ 
1 -2v 



[3M iso c - V • u] . 



(17) 



(18) 



If we take the curl of Equation (|17p . we obtain V 2 (V x u) = 0. From potential theory, we know 
that u must then be curl free and can be written in terms of the gradient of a scalar potential: 
u = V^. This gives the equation 



V 



1 + v 
l-v 



0. 



(19) 



This can then be integrated with the constant of integration taken as zero without loss of generality. 
We hence have the dilatation, V • u, as 



M iso ( I r 



l — u / 

and can use this to rescale the enthalpic contribution in the chemical potential as 



li = a — 2bc + ksT log 



(20) 



(21) 



Note that this result was first obtained and analyzed for phase-field models in [4iy42j. 

We can next obtain a spinodal region in the parameter space, which is the region in which 
homogeneous concentrations are linearly unstable. The homogeneous state of the system can be 
calculated by taking c = cq to be constant, and in perturbing this state with the infinitesimal 
quantity c = c — cq , the linearized evolution equation for the concentration perturbation will take 
the form 



dc 
dt 



c V • <^ BV 



-26 + 



c (l - c ) 



c - V • (KVc) 



(22) 
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Figure 1: Spinodal region in terms of the isotropic lattice mismatch M- mo and homogeneous con- 
centration Co- The stable and unstable regions are labeled within the figure. The parameters used 
are a = 1.6 x 10" 20 J 03], p = 14nm- 3 , G = 5x 10 10 Pa, v = 0.22 [33], and k B T = 4 x 1(T 21 J at 
room temperature. 



By assuming the normal mode of the perturbation as c = cexp(at + iq • x) 
of a given wave vector q, this gives the dispersion relation 



cr(q) = - Co q • (Bq) 



-26 + 



k B T 



co(l - c ) 



q • (Kq) 



where a is growth rate 



(23) 



The spinodal can now be calculated for parameters in which there exists a positive growth rate, 
and hence we have the region defined by 



M 2 < 



1 



1 



2a 



k B T 



co(l - c ) 



(24) 



It is important to note that this stability criterion is independent of the sign of Mi so , and hence both 
crystal expansion and contraction will have the same effect. Additionally, the increase of elastic 
strains narrows the stability region in parameter space and can even annihilate the spinodal for a 
given interaction parameter a as shown in Figure (pQ). The parameters obtained from experiment 
are listed in the caption. This latter point is somewhat troubling as not only does Li^FePC^ exhibit 
both high strains and phase-separation upon intercalation, but the phase-separation in this material 
has been observed to occur along the direction of greatest strain [15|I44| . This apparent contradiction 
is resolved in the following section. 



4 Anisotropic Coherency Strain 

To account for the alignment of the FeP04/LiFeP04 phase boundary with the (b, c) plane of the 
host crystal, we must include anisotropic properties of the system into our model (note that we 
are using the standard notation for the axis labels of Li x FeP04 as seen in [15]). In particular, a 
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more appropriate approximation for the mismatch tensor corresponding to this system would be 
diagonal with distinct values. For L^FePC^, these diagonal elements are roughly: Mn = 0.052, 
M22 = 0.036 and M33 = —0.019 [15J. To capture the effects of anisotropic coherency strains, we 
will simplify the model by reducing the remaining tensors to scalars as B = BI and K = KI, as 
they merely act as spatial scalings. We can then easily nondimensionalize the spatial and temporal 
coordinates with the following relations 



A 2 x 

x — > Ax, t — > —t, A 
' V 



K 
kef' 



(25) 



where we have used the Einstein mobility relation to introduce the diffusivity T> = ksTB. The 
characteristic length-scale, A, represents the typical width of the FeP04/LiFeP04 interface which 
is on the order of nanometers 1151. We hence obtain the concentration evolution 




-2ac + log 



1 -c 

3 3 



3 3 



i=l j=l 



V 2 c--l—YYM ll C ij ^- 



where a 



(26) 
(27) 



which is coupled with the equilibrium conditions (jlOp . To further further illuminate the effects of 
coherency strain anisotropy, we again take the host material to be an isotropic body. While this is 
not quite physical, simulations have been performed for both isotropic and anisotropic bodies using 
the elastic moduli found in |33| to verify this approximation, and no qualitative differences have 
been observed. Furthermore, as also seen in [33], the anisotropy is much stronger in the lattice 
mismatch (e.g. a contraction along the c-axis is induced by insertion), and hence these effects are 
amplified by the fact that while the energetic contributions from elasticity are (D(Cij), they are 
indeed quadratic in the mismatch parameters [35]. This system therefore simplifies to 

^ = V • jcV -2ac + log f j-^J - v2c - 2 7 ( MV + 7 _ ~^7; tr { M } V ) ' u } ' ( 28 ) 
(1 - 2u) V 2 u + V(V • u) = 2 [(1 - 2^)MV + v tr{M}V] c, (29) 



where a 
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tr{M 2 } + 



1 - 2v 



tr{M} s 



G 



7 



pk B T 



(30) 



and tr{o} denotes the trace operator. As before, we perturb our system about the basic state of 
constant concentration and zero deformations to obtain the linearized system 



dc 



cqV 2 



-2a 



1 



co(l - c ) 



c - V 2 c - 2 7 I MV + 



1 - 2v 



tr{M}V • u 



(1 - 2v) V 2 i + V(V • u) = 2 [(1 - 2i^)MV + v tr{M}V] c. 



(31) 
(32) 



By introducing the vector of fields v = (ui, U2, U3, c), we can again take the normal mode of the 
perturbations as v = vexp(at + iq • x) and obtain a linear algebraic system of the form Av = 0. 
This operator can then be written as A(<r, q) = 5](cr) + Q(q), where £44 = a is the only nonzero 
element of X, and the elements of Q are given in the Appendix. For the perturbation eigenvector 
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v to have nontrivial solutions, the operator A must be singular (i.e. det{A} = 0), and hence we 
obtain the dispersion relation for the growth rate of a given spatial mode 



-det{Q(q)} 
det{Q(4)(q)} 



(33) 



where the parenthetical superscript denotes the submatrix of Q obtained by removing the fourth 
row and column. A typical plot of a can be seen in Figure ([2]). For presentation purposes, we have 
only shown a slice in the (x, y) plane. Lighter colors denote larger growth rates, the curve represents 
the zero (neutral) contour of a, and the x's mark the points of greatest growth. The parameters 
are the same as those used in Figure (pQ) with the addition of lattice mismatch rates listed in the 
caption obtained by [15J. Note that the most excited wave- vectors lie along the x-axis, which for 
the parameters used, is the direction of greatest strain. To further understand this, we take the 
simple case of unidirectional strain {M%\ = Af U ni> M22 = -^33 = 0). If we then write the growth 
rate as <r(q) = oo(q) + °"s(q)> where gq is the growth in the absence of strain (i.e. M un ; — > 0), and 
a s is the strain induced contribution, then this contribution takes the form 



It can be seen here that this contribution is minimized for a given wave vector of magnitude 
q in the {q y ,q z ) plane, which means that phase-separation is suppressed the most in directions 
orthogonal to the unidirectional strain. Additionally, along the vector q = (q x ,0,0), the strain 
induced contribution is maximized at a s = 0, which means that the direction of strain for this 
case is the only direction without suppression. It can be shown that if all Ma have the same sign 
(i.e. pure expansions or pure contractions), then phase-separation will always occur in the direction 
in which \Ma\ > \Mjj\, {M^kl (i.z. the direction of largest strain), however it will be seen in the 
next section that the presence both contractions and expansions result in the possibility of skewed 
phase interfaces. 

5 Numerical Simulations 

To verify the results from the linear stability analysis, we perform numerical simulations of the 
governing equations (l28M29p . as this analysis can only predict the initial dynamics still within the 
linear regime. Both two- and three-dimensional simulations are presented to represent two limits 
of this system. 

5.1 Phase Separation with Anisotropic Plane Strain 

We first present the limit in which strains along the c-axis are neglected (i.e. U3 — > 0) under 
the assumption that Li transport along molecular channels direction is in constant equilibrium 
due to the very fast mobility p2] and suppressed phase separation [23]. This is a reasonable 
assumption for reaction-limited Li x FeP04nano particles [19j|22], but would break down in larger 
particles (> lOOnm) due to channel-blocking Li/Fe anti-site defects, which lead to size-dependent 
effective mobility [TJ]. Note that the stability results from section @ remain unchanged for the 
two-dimensional case. For numerical convenience, we will further approximate the flux (before the 
rescaling) as j = —BVfi. This is a common assumption in phase-field modeling since the prefactor 
c merely scales the initial time evolution and does not affect long-term behavior of the system which 




q = |q|- 



(34) 
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-1 1 



Figure 2: Dispersion relation with anisotropic mismatch strains. The curve indicates neutral sta- 
bility (i.e. a = 0), while the x's mark the largest growth rates. The lattice mismatch parameters 
used are M n = 0.052, M 22 = 0.036 and M 33 = -0.019 [33]. 

is solely determined by the form of the chemical potential [30,35,41]. We use a pseudo-spectral 
method for the spatial differentiation with 256 x 256 modes and a semi-implicit scheme for temporal 
updating. The linear terms are handled with the backward Euler method, which although only first 
order accurate in time, provides the necessary stability required for the stiff solutions of the system. 
The second order Adams-Bashforth method is used for the nonlinearities, and periodic boundary 
conditions are used to approximate the infinite domain. As the equations for the deformation fields 
are linear, they can be solved exactly in Fourier space in terms of the order parameter, and hence 
the system can be reduced to a single evolution equation. If we define the Fourier transform of the 
order parameter as J{c(x, t)} = 0(q, t), where q is the corresponding wave-vector, the evolution of 
(f) in Fourier space will be therefore governed by 

£ - - ' 4 + * - rt {* (t^) } • (35) 

where the contribution from the coherency strains is given as 

M(q) = [(1 - v)M n + vM 22 \ [((1 - v)q 2 + q 2 ) M u + (vq 2 - q 2 ) M 22 ] q 2 x (36) 
+ [pM n + (1 - v)M 22 \ [{vq 2 - ql) M n + ((1 - u)q 2 + q 2 x ) M 22 ] q 2 y . (37) 

A typical simulation evolved from the state c = 0.5 perturbed with small amplitude noise is 
shown at various times in Figure ((3j). These perturbations are initially damped, and regions of 
localized concentrations begin to form. The regions coarsen until the concentration is aligned with 
the y-axis confirming the spontaneous phase-separation along the direction of greatest strain. Note 
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t = t = 1 




02468 02468 

x x 



Figure 3: Time evolution of system (|28H29f) starting from the state c = 0.5 perturbed with small 
amplitude noise at t = 0, 1, 7, 10 as labeled. 

that the parameters are the same as those listed in the captions of Figures ([1]) and ([2]) excluding 
contributions from the third dimension. It should be mentioned that this model is quite robust, as 
similar results can be obtained from simulations over a wide range of parameters. 

5.2 Three Dimensional Coherency Strain in Li x FePO*4 

For a more complete picture of phase separation in Li x FeP04, it is necessary to consider fully 
anisotropic three-dimensional coherency strain. In particular, we perform full three-dimensional 
simulations to include the effects of contractions in the c-direction. As seen in Figure Q, the 
equilibrium state results in skewed interfaces, which have also been seen in some experiments |lipi3j . 

6 Conclusions 

We have considered the intercalation of Li within Li x FeP04 as a model system which exhibits both 
phase-separation and highly anisotropic coherency strains. A phase-field model which incorporates 
energetic contributions from entropy, enthalpy and elastic properties of the host material has been 
developed. We have shown through linear stability analysis that while coherency strains suppress 
phase-separation, this suppression is maximized in directions orthogonal to the direction of a given 
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Figure 4: Equilibrium state of Li^FePC^with anisotropic coherency strain in three dimensions 
showing skewed, striped phases. 

strain. Hence the system is expected to phase-separate along the direction of greatest strain (a- 
axis) given sufficiently large enthalpic contributions. As linear stability theory can only predict the 
initial evolution of spinodal decomposition, numerical simulations were performed on an idealized 
system to illuminate the effects of coherency strain anisotropy, and phase boundaries were observed 
to align perpendicularly to the a-axis as observed in experiment. Coherency strain anisotropy has 
therefore been shown to be a possible mechanism for this alignment. We have also considered 
phase separation in the three dimensional case with fully anisotropic coherency strain, including 
contraction along the c axis. This leads to skewed phase boundaries, which have also been seen 
in some experiments. These interesting observations are unified and extended to electrochemically 
driven intercalation at fixed current in Ref. |23j. 

Acknowledgements 

This work was supported by the National Science Foundation under Contracts DMS-0842504 and 
DMS-0948071 and partially under the auspices of the US Department of Energy by LLNL under 
Contract DE-AC52-07NA27344. 



11 



Appendix 

The elements of Q are listed as follows, where we have defined q = |q| with q = (q x ,Qy,Qz)- 

Qn = (1 - 2v)q 2 + q 2 x , Qi2 = q x q y , Qi3 = QxQz, Qi4 = 2i[(l - 2v)M u + vM kk \q x , 
Q2i = q y q x , Q22 = (1 - 2v)q 2 + q 2 , Q 23 = q y q z , Q24 = 2i[(l - 2u)M 22 + vM kk ]q y , 
Q3i = q z qx, Q32 = q z q y , Q33 = (1 - 2v)q 2 + q 2 z , Q u = 2i[(l - 2u)M 33 + vM kk ]q z , 
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